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ABSTRACT 

Data taken in stereo mode by the High Resolution Fly's Eye (HiRes) air fluorescence experiment 
are analyzed to search for correlations between the arrival directions of ultra-high-energy cosmic rays 
with the positions of BL Lacertae objects. Several previous claims of significant correlations between 
BL Lacs and cosmic rays observed by other experiments are tested. These claims are not supported 
by the HiRes data. However, we verify a recent analysis of correlations between HiRes events and a 
subset of confirmed BL Lacs from the 10th Veron Catalog, and we study this correlation in detail. 
Due to the a posteriori nature of the search, the significance level cannot be reliably estimated and 
the correlation must be tested independently before any claim can be made. We identify the precise 
hypotheses that will be tested with statistically independent data. 
Subject headings: cosmic rays — BL Lacertae objects: general — galaxies: active 




I. INTRODUCTION 

Among the most striking astrophysical phenomena to- 
day is the existence of cosmic ray particles with energies 
up to and exceeding 10^" eV. It is currently unknown 
where and how these particles are accelerated to such en- 
ergies and how they travel astronomical distances with- 
out substantial energy loss. 

In an attempt to understand the origin of these 
particles, the limited world data set of ultra-high- 
energy cosmic ray arrival directions has been sub- 
jected to extensive searches for correlations with the 
positions of objects from known astrophysical source 
classes. In particular, significant correlations be- 
tween subsets of BL Lacertae objects and cosmic ray 
arrival directions recorded by the Akeno Giant Air 
Shower Array (AGAS A) and the Yakutsk expe r iment 
have been claimed fTinvakov fc TkachevI I2001L 120021 
iGorbunov et alJ I2002D . Searches for correlations with 
BL Lac objects are well-motivated. BL Lacs are a sub- 
class of blazars, which are active galaxies in which the 
jet axis happens to point almost directly along the line 
of sight. The EGRET instrument on board the Comp- 
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ton Gamma Ray Observatory (CGRO) has firmly es- 
tablished blazars as sources of high energy 7-rays above 
100 MeV l|Hartman et al.lll999j) . and several BL Lac ob- 
jects have been observed at TeV ener gies with ground- 
based air Cherenkov telescopes f see iHoran fc Wccke^ 
12004 for a summary). High energy 7-rays could be by- 
products of electromagnetic cascades from energy losses 
associated with the acceleration of ultra-high-energy cos- 
mic rays and the ir propagation in intergalactic space 
( Bcrczinsk ii et al.l ir990: Covvi & Aharonian 1997). 

However, the claimed correlations between ultra-high- 
energy cosmic ray arrival directions and BL Lac ob- 
jects are controversial. A problematic aspect of the 
claims is the procedure used to establish correlations and 
evaluate their statistical si gnific ance. Several authors 
l|Tinvakov fc Tkachevll20?)lIl20(illGorbimov et alll200^ 
explicitly tuned their selection criteria to assemble cat- 
alogs that show a maximum correlation with arrival di- 
rections of cosmic rays above some energy. An unbiased 
chance probability for these correlations can then only 
be arrived at if the claim is tested on a statistically in- 
dependent data set. Since the available data set is small, 
this rigorous procedure is often abandoned, and instead 
an attempt is made to estimate a statistical penalty fac- 
tor to compensate for the number of trials involved in 
the tuning. "Hidden" trials, unfortunately, make these a 
posteriori estimates highly unreliable, and claims of BL 
Lac correlations have been criticized on th ese grounds 
l|Evans et al.ll20mt ISteT^ fc PoutanenI 120051) . Addition- 
ally, in some cases it has been shown that statistically 
independent data sets do not confirm the correlations 
(Torres et al. 2003). 

The operation of the stereoscopic High Resolution 
Fly's Eye (HiRes) air fluorescence detector is providing 
a large data set of cosmic ray events with unprecedented 



2 



The HiRes Collaboration 



angular resolution for the study of small-scale anisotropy 
and source correlations. In this paper, we report on 
searches for correlations between BL Lac objects and 
HiRes stereoscopic events observed between 1999 Decem- 
ber and 2004 January. The quality cuts applied to this 
data sample are described in detail in lAbbasi et alJ l)2004 
|2005). 

The outline of the paper is as follows: in Section |21 we 
describe our use of the maximum likelihood method for 
correlation searches with multiple sources. In Section |3| 
we consider the previous claims of correlations between 
BL Lacs and ultra-high-energy cosmic rays and test them 
using HiRes stereo data. In Section 21 we perform a 
general search for correlations between the HiRes data 
set and confirmed BL Lacs. In Section [SJ we summarize 
the results and discuss further studies. 

2. MAXIMUM LIKELIHOOD METHOD 
2.1. Description 

To address a number of shortcomings of binned anal- 
yses of cosmic ray arrival directions, we have recently 
applied an unbinned maximum likelihood method in the 
search for point sou rces of ultra-high-energy cosmic rays 
dAbb asi et al.ll2005(l . This approach uses the probability 
density function for each individual event rather than 
requiring a fixed bin size. Two important advantages 
of this method are the ability to accommodate events 
with different errors, and to give weighted sensitivity 
to angular separations — avoiding the loss of information 
that follows from choosing an angular separation cut-off. 
With minor modifications, the same maximum likelihood 
method can also be used to search for correlations with 
a specified hst of potential sources. 

The premise involved in the maximum likelihood anal- 
ysis is that the data sample of N events consists of 
source events which came from some source position(s) in 
the sky, and N — rig background events. A background 
event arrives according to the probability distribution 
given by the detector exposure to the sky, i?(x), where 
X are equatorial coordinates. The true arrival direction 
of a source event is the location of the source s, but 
the event is observed somewhere near s according to the 
probability distribution Qi(x, s), where Qi depends on 
the angular uncertainty in arrival direction of the ith 
event. 

Because it is not known whether a given event is a 
source or background event, the probability distribution 
function (or "partial" probability) for the ith event is a 
weighted sum of the source and background probability 
distributions: 

T) J\ — Ti 

P,(x) = ^Q,(x,s) + ^^i?(x) . (1) 

This describes the distribution of arrival directions un- 
der a single-source hypothesis. For a hypothesis with M 
sources, we must modify Q to include multiple source lo- 
cations. We will assume for this analysis that the sources 
have equal luminosity. In this case, we only need to com- 
pensate for the varying exposure of the detector to dif- 
ferent parts of the sky: the probability for a source event 
to come from the jth source is proportional to the detec- 
tor exposure R{sj) to the source location Sj. The total 
source probability distribution for the ith event is 



then the weighted sum of the individual source probabil- 
ities: 

M M 

Qr(x)-5]i?(s,)g.(x,s,)/5]i?(sfe) . (2) 

j=i k=i 

Replacing Qi in Eq.^with we evaluate the partial 
probability of the ith event at its observed location x^: 

P,(x,) = ^gr(x.) + ^^^i?(x.) . (3) 

The likelihood for the entire set of N events is then 
given by: 

N 

^ns)=l[P^{^^) . (4) 

i=l 

The best estimate for the number of events contributed 
by the sources is determined by finding the value of 
that maximizes C. 

Because the value of the likelihood function depends 
on the number of events, a more useful quantity than C 
is the likelihood ratio TZ: 

where C{0) is the likelihood function of the null hypoth- 
esis (tIs = 0). In practice, we maximize In 7?., which is 
equivalent to maximizing C. 

The significance of the resulting In TZ can be estimated 
using — 2 In??.. When is positive, this agrees well 
with the distribution for one degree of freedom. Be- 
cause Us corresponds to the excess number of events cor- 
relating with source positions, a negative best-fit value 
for Us will occur whenever there are fewer events near 
source positions than expected. Negative values are 
not physically meaningful in the point-source search, but 
they are useful for evaluating significances. To distin- 
guish an excess in correlations from a deficit, we assign 
the negative solution x = — \/2 In 7?. when the best-fit 
is negative. We can then check the significance estimated 
from the x^ distribution by performing the same likeli- 
hood analysis on simulated data sets and ranking them 
according to their x values. We will use T to denote the 
fraction of simulated, isotropic event sets which yield a 
value of X greater than or equal to that of the data. 

2.2. Implementation 

For the source probability function Qi we employ a 
circular Gaussian of width ai corresponding to the an- 
gular uncertainty of the ith event, as estimated by the 
stereo event reconstruction. The mean of the angular un- 
certainty of the HiRes stereo events is slightly larger at 
lower energies, growing from (cr) = 0.44° above 10^* eV 
to (cr) = 0.60° below lO" '^^ eV. This of course is ac- 
counted for by the use of individual errors in the maxi- 
mum likelihood analysis. 

Because the geometrical acceptance of the detector is 
a function of energy, we determine the background prob- 
ability distribution i?(x) in two different ways. For large 
event samples (> 1000), we generate a background dis- 
tribution from a full time-swapping of the data itself: 
the equatorial coordinates of each event are recalculated 
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TABLE 1 

Previously claimed correlations between BL Lacs and cosmic rays, and tests with independent HiRes data 



Sample 


CR Data and Energy 


# 


Binned Analysis 


Max. Like. Analysis 


(# Objects) Reference 


Threshold (EeV) 


Ev. 


Bin Size Pairs Prob. 


InTZ ris 


T 



Claim 1: 


A (22) 


TTOl 


AGASA>48, Yak.>24 


65 


2.5° 


8 


< 10-4 








Test: 




HiRes > 24 


66 


2.5° 





1.00 


(0) 


(0) 


0.75 


Claim 2: 


B (157) 


TT02 


AGASA > 40 


57 


2.5° 


12 


0.02 








Test: 




HiRes > 40 


27 


2.5° 


2 


0.78 


(0) 


(0) 


0.26 


Claim 3: 


C (14) 


G02 


AGASA>48, Yak.>24 


65 


2.9° 


8 


10-4 








Test: 




HiRes > 24 


66 


2.9° 


1 


0.70 


(0) 


(0) 


0.68 



References. — [TT Ol] jTinvakov fc TkachevI 1200111 : [TT02] ITinvakov fc Tkache-J 1200211 : [G02] IGorbunov et al.l 1200211 : [G04] 
Ig^ unov et alj 120041 . In [TT02] and [G02], the authors also attempt to correct for the deflections of charged primaries by the 
galactic magnetic field; these results are not considered here. 



using all of the recorded event times, and i?(x) is the 
cumulative map of all of these virtual event locations 
convolved with a circular Gaussian function for smooth- 
ing. For small event samples (e.g. the 271 events above 
10^^ eV) the data set is too sparse to generate a useful 
time-swapped map. Instead, we rely on a full detector 
simulation to estimate the local geometrical acceptance, 
and convolve this with the event times to generate i?(x). 
This procedure a nd the detector simu lation are described 
in more detail in lAbbasT "eTall l|200|) . 

3. TESTS OF PREVIOUS CORRELATIONS 
OBSERVED WITH AGASA AND YAKUTSK 
DATA 

We briefly review some past claims of correlations be- 
tween cosmic ray arrival directions and BL Lacs, and 
then test these with HiRes data under the same con- 
ditions. All samples of BL Lacs are selected from 
objects classified as "BL" in Table 2 of the Veron 
Catalog of Quasars and AGN, 9th or 10th Editions 
ijVeron-Cettv Veroull2?M l2nnT|> . 

• Samp le A: described in ITinvakov fc TkachevI 
l|20nH ). contains 22 BL Lacs from the Veron 9th 
Catalog selected on the basis of optical magnitude 
(to < 18), redshift (z > 0.1 or unknown), and 6cm 
radio flux (Fg > 0.17Jy). 

• Samp le B: described in ITinvakov fc TkachevI 
l)2002f) . contains 157 BL Lacs from the Veron 10th 
Catalog with optical magnitude m < 18. 

• Sample C: described in IGorbunov et alJ l)2002j) . 

consists of 14 BL Lacs from the Veron 10th cat- 
alog that were selected by the authors on the basis 
of possible association with identified and uniden- 
tified 7-ray sources in the Third EGRET Catalog 
llHa,rtma,n et a].llT99l. 

Table ^ shows the correlations originally claimed us- 
ing these BL Lac samples and cosmic ray data from the 
AGASA and Yakutsk experiments. The energy thresh- 
olds and angular bin sizes vary from analysis to analysis 
as shown. The results of testing each claim as nearly as 
possible with an equivalent set of HiRes data are also 
presented. Both a binned analysis with the originally 
used bin size and a maximum likelihood analysis using 



the point spread function of individual HiRes events are 
performed. In the binned analysis, the number of event- 
object pairs with angular separation less than the bin size 
are counted, and the probability for the same or greater 
number of pairs is evaluated using simulated isotropic 
event sets. None of the three previous claims of correla- 
tions based on other cosmic ray data sets are confirmed 
by the tests. Each test, in fact, finds a deficit or no excess 
of HiRes events correlating with BL Lacs, indicated by 
(0) values for In TZ and rig . The fraction !F of simulated 
sets with stronger correlation than the data is calculated 
as described above. 

In the tests of Claims 1 and 3, the size of the HiRes 
event sample is comparable to the size of the combined 
AGASA and Yakutsk event samples. Assuming a Poisson 
distribution with mean number of event - BL Lac pairs 
given by the AGASA and Yakutsk results, the observed 
number of HiRes - BL Lac pairs excludes the claimed 
correlations at a confidence level greater than 99 % in 
each case. In the test of Claim 2, the HiRes event sam- 
ple is smaller than that of AGASA. Here the claimed 
correlation is excluded at the 90 % confidence level. 

4. RECENT CORRELATIONS OBSERVED WITH 
HIRES DATA 

Recently, the p ublished HiRes e vents above 10"'^^ eV 
were analyzed bv IGorbunov et all (j2004), and correla- 
tions with the BL Lacs of Sample B were claimed at 
the 10"'^ level. The analysis used a fixed bin size of 
0.8°, which the authors argued is optimal for a point- 
source search given the HiRes angular resolution. We 
verify this analysis by applying the maximum likelihood 
method to the same data set and source sample, and find 
InTZ — 6.08 for Ug — 8.0; the fraction of Monte Carlo sets 
withjugherjii^js = 2 x 10"^. 

IGorbunov et al.l 1)2004(1 analyzed the entire set of 
HiRes events above 10^^ eV because the individual 
event energies were not published. Therefore this en- 
ergy threshold was not tuned to maximize correla- 
tions with BL Lacs. However, because the original 
claim (,Tinvakov fc Tkachev-2002i) was based on AGASA 
da ta with energies above 4 x 10^^ eV, the correlation 
in IGorbunov et al.l l)2004|) does not confirm a previous 
claim, but rather represents a new hypothesis. This is 
demonstrated by the fact that the HiRes data shows no 
excess correlation with this sample of BL Lacs when the 
same 4 x 10^^ eV energy threshold is used, as indicated 
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18.6 18.8 19 19.2 19.4 19.6 19.8 20 
log(Energy/eV) 

Fig. 1 . — In TZ result as a function of minimum energy threshold of 
the HiRes data set. The 10^^ eV energy threshold of the published 
data is indicated. 



in the test of Claim 2 in Tabled 

The observed correlation warrants further investiga- 
tion. We report on extending the analysis to lower en- 
ergy HiRes data and to the rest of the confirmed BL Lacs 
in the Veron catalog. 

4.1. Event Sample: Energy Dependence of Correlations 

An important question is whether and how the ob- 
served correlation depends on the energy threshold. Fig- 
ure ^ shows the result of the same maximum likelihood 
analysis above, performed repeatedly using increasing 
energy thresholds from lO^^-^ eV to lO^" eV. The 10^^ eV 
threshold corresponding to the published data set is in- 
dicated, and it clearly stands out as the threshold that 
gives a local maximum in the significance of the correla- 
tion. 

One of the motivations for using an energy threshold 
in small-scale anisotropy searches is that charged cosmic 
ray primaries are subject to deflection by galactic and 
extragalactic magnetic fields. The highest energy pri- 
maries will be least deflected, and consequently will be 
the best candidates for correlation studies. However, an 
a priori energy threshold cannot be established, because 
detailed knowledge of the galactic and extragalactic mag- 
netic fields is lacking. 

Figure H indicates that most of the correlation comes 
from events with energies between 10^^ eV and 10^^'^ eV. 
At these energies, it is generally assumed that the Galac- 
tic magnetic field will deflect a proton primary by many 
degrees; nuclei will be deflected even more. In spite 
of this, the correlations are consistent with the ~ 0.5° 
scale of the detector angular resolution. This might im- 
ply that the correlated primary cosmic rays are neutral, 
thus removing the motivation for restricting the analysis 
to events with energies above some arbitrary threshold. 
A search for correlations with the entire HiRes stereo 
data set of 4495 events at all energies is justified. 

Applying the analysis to the entire data set and Sample 
B gives ln7^ = 6.16 for = 31, with T = 2x lO'^. 
This of course includes the effect of the correlated events 
above 10^^ eV; for the independent sample of 4224 events 
below lO^^eV, we find ln7^ = 3.10 for 22, with 

jF = 6x 10-3. 

4.2. Source Sample 

The Sample B of BL Lacs discussed above consists of 
157 confirmed BL Lacs in the 10th Veron Catalog with 
optical magnitude m < 18 and which are classified as 
"BL" . The rest of the confirmed BL Lacs are classified 
as "HP" (high polarization). It is natural to ask about 



TABLE 2 

HiRes — BL Lac Correlation Results: Fraction 

:r OF SIMULATED HiReS SETS WITH STRONGER 
CORRELATION SIGNAL. 



Source Sample {# Obj.) 


All Energies 


E >10 EeV 


"BL" (157) 


2 X 10"" 


2 X 10"" 


"HP" (47) 


0.3 


6 X 10-3 


"BL"-|-"HP" (204) 


5 X 10""' 


10-5 



Note. — Correlations are with confirmed BL Lacs in 
Table 2 of the Vcron 10th Catalog IVeron-Cettv &: VeronI 
1200111 . classified as cither "BL" or "HP," with m < 18. 



these objects as well. Indeed, of the six blazars which 
have confirmed detections in 7-rays at TeV energies, half 
are classified as "HP" and half as "BL" in the Veron 
catalog. 

We apply the same cut on optical magnitude m < 18 as 
in Sample B to the "HP" objects, and arrive at a sample 
of 47 objects. The result of the maximum likelihood 
analysis applied to this independent sample of BL Lacs 
and the HiRes events above 10^® eV is \nTZ = 3.13 for 
Us = 3.0, with = 6 X IQ-^. We also perform the same 
analysis on the events below 10^^ eV. No excess is found. 

The results for HiRes events of energy above 10^^ eV 
and all energies with BL Lacs classified as "BL," "HP," 
or "BL" and "HP" combined are summarized in Table |21 

The equivalent analyses have been performed on the 
same classes of BL Lacs with m > 18: no excess corre- 
lation is found in any of these cases. It is apparent from 
th ese results that the m < 18 cut which was identified 
in iTinvakov fc TkachevI l|2002D as optimal for AGASA 
also isolates the BL Lac objects which show excess cor- 
relations with HiRes events. Under the BL Lac source 
hypothesis, of course, it is not unreasonable to expect 
the closer and more luminous objects to contribute more 
strongly. However, since the Veron catalog is not a uni- 
form sample of BL Lac objects, the interpretation of this 
cut may involve a more complicated interplay of selection 
effects from the underlying surveys which make up the 
catalog.^ 

4.3. TeV Blazars 

Among the closest and brightest of the "BL" and "HP" 
BL L acs are six whi ch are confirmed sources of TeV 7- 
rays ijHoran fc Weekes .2004) . Five of these, shown in 
Table |31 are high in the northern sky and well within 
the field of view of HiRes. We perform the maximum 
likelihood analysis on this set of objects using all of the 
HiRes data, and find \nTZ = 4.78 for Ug — 5.6 with — 
1Q~^. For just the HiRes events above 10^^ eV, the result 
is ln7^ = 6.15 for = 2.0 with = 2x10-^. In TableEl 
we show the results for treating each blazar in turn as a 
single source hypothesis. 

^ The Veron catalog strives to be "complete" only in the sense 
of a complete survey of the literature and catalog of all known BL 
Lacs; it does no t represent an unbiased sta tistical sample of BL 
Lacs in any way IVeron-Cettv fc Veronl2000l I2OOII) . This does not 
exclude the possibility of using subsets of the catalog to identify 
correlations with cosmic rays, but it means that any inferences 
about the BL Lacs based on such correlations may be highly biased 
and simply an artifact of the underlying combination of different 
surveys in the catalog. 
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TABLE 3 

TeV Blazar Correlation Results with HiR.es events 
(all energies) 



TeV Blazars 




Correlation Results 


Name 


z 


V Mag 


ris 


InTl 


T 


Mrk 421 


0.03 


12.9 


0.3 


0.04 


0.2 


H1426+428 


0.13 


16.5 


(0)- 


(0) 


0.4 


Mrk 501 


0.03 


13.8 


3.3 


5.27 


6 X lO""' 


1ES1959+650 


0.05 


12.8 


2.0 


2.87 


8 X 10-3 


1ES2344+514 


0.04 


15.5 


{0)b 


(0) 


0.7 


Combined Set 






5.6 


4.78 


10^3 



^No excess: < 3.5 at 90% confidence level. 
'^No excess: < 2.4 at 90% confidence level. 



TABLE 4 

HiRes — BL Lao Correlation Summary: Fraction T of 
simulated HiRes sets with stronger correlation signal. 



Source Sample (# Obj.) 


All Energies 


E > 10 EeV 


"BL" Objects, m < 18 (157) 


2 X 10-"* 


2 X 10-* 


Confirmed BL Lacs, m < 18 (204) 


5 X 10-"' 


10-5 


Confirmed TeV Blazars (6) 


10-3 


2 X 10-* 



Note. — All samples arc contained within Table 2 of the Vcron 10th 
Catalog. The samples overlap and arc not independent: "Confirmed 
BL Lacs" combines "BL" and "HP" classified BL Lacs; TeV Blazars 
are a subset of the confirmed BL Lacs. 



5. RESULTS AND DISCUSSION 

In this paper, we have used an unbinned maximum 
hkchhood method to analyze correlations of ultra-high- 
energy cosmic ray arrival directions with BL Lac ob- 
jects in the Veron 10th Catalog. We first tested previous 
claims of correlations between BL Lacs and cosmic rays 
which were based on AGASA and Yakutsk data. Using 
the independent HiRes stereo data set, these correlation 
claims are excluded at the 99 % (Claims 1 and 3 in Ta- 
ble or 90% confidence level (Claim 2). 

H owever, we h ave verified the observation 
by iCorbunov et all l|2004D that the set of HiRes 
stereo events with energies above 10^^ eV shows an ex- 
cess of events correlated with confirmed BL Lacs marked 
as "BL" in the Veron 10th Catalog. We emphasize that 
the observed correlation does not confirm a previous 
claim, because it requires a lower energy threshold. It 
can only be confirmed with new data. 

We have explored the extension of the analysis to 1) 
HiRes events of all energies, and 2) the rest of the con- 
firmed BL Lacs (labeled "HP") in the Veron 10th Cata- 
log. In each case, correlations at the significance level of 



^ 0.5% are found. While statistically independent from 
the above result, these are not strictly tests of that claim. 
However, the combination offers well-defined hypotheses 
which can be tested with new data. 

The results of combining the analysis of low and high 
energy events and "BL" and "HP" BL Lacs are summa- 
rized in Table 01 Also shown are the results for HiRes 
events and the subset of BL Lacs which are confirmed 
sources of TeV 7-rays. 

The analyses described here have only been performed 
on the data recorded through 2004 January. The HiRes 
detector will continue observations through the end of 
2006 March. By that time the independent sample of 
data since 2004 January is expected to reach approxi- 
mately 70 % of the size of the sample analyzed here. This 
will provide an opportunity to test the correlations in Ta- 
ble^ We note that while the correlation signals appear 
stronger for the events above 10^^ eV, a conservative ap- 
proach which includes consideration of the entire data set 
will avoid the possibility that a real correlation has been 
"over-tuned" by an arbitrary threshold and is missed in 
a future analysis. 

As mentioned earlier, real correlations on the scale of 
the detector angular resolution would suggest neutral 
cosmic ray primaries for these events, or at least that 
the primaries were neutral during significant portions 
of their journey through galactic and extragalactic 
magnetic fields. Primaries such as neutrons and photons 
are problematic, however, because of short mean free 
paths (~ a few Mpc) at these energies. As a fluorescence 
detector, HiRes can study the height of the maximum 
of the shower development, a parameter that has 
some sensitivity to the primary particle type. Showers 
induced by photons, for example, tend to develop lower 
in the atmosphere than those induced by nucleons. An 
analysis is underway to address the question whether 
the correlated HiRes stereo events are compatible with 
gamma-induced showers. 
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